Bond algebraic liquid phase in strongly correlated multiflavor cold atom systems 
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When cold atoms are trapped in a square or cubic optical lattice, it should be possible to pump 
the atoms into excited p— level orbitals within each well. Following earlier work, we explore the 
metastable equilibrium that can be established before the atoms decay into the s— wave orbital 
ground state. We will discuss the situation with integer number of bosons on every site, and consider 
the strong correlation "insulating" regime. By employing a spin-wave analysis together with a new 
duality transformation, we establish the existence and stability of a novel gapless "critical phase", 
which we refer to as a "bond algebraic liquid" . The gapless nature of this phase is stabilized due to 
the emergence of symmetries which lead to a quasi-one dimensional behavior. Within the algebraic 
liquid phase, both bond operators and particle flavor occupation number operators have correlations 
which decay algebraically in space and time. Upon varying parameters, the algebraic bond liquid 
can be unstable to either a Mott insulator phase which spontaneously breaks lattice symmetries, 
or a Z2 phase. The possibility of detecting the algebraic liquid phase in cold atom experiments is 
addressed. Although the momentum distribution function is insufficient to distinguish the algebraic 
bond liquid from other phases, the density correlation function can in principle be used to detect 
this new phase of matter. 

PACS numbers: 75.10.Jm, 71.10.Hf, 05.30.Jp, 03.75.Lm 
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I. INTRODUCTION 

Interacting bosonic systems hopping on two or three 
dimensional lattices (d = 2, 3) typically have one of two 
ground states: superfluid or a Mott insulator—. Search- 
ing for bosonic phases other than these two has attracted 
much attention, especially gapless bosonic phases which 
do not break the global U(X) symmetry. In such gap- 
less featureless bosonic phases, the correlation functions 
between physical operators usually fall off algebraically. 
Being the lattice boson analog to the algebraic spin liq- 
uids which can occur in spin models, these phases can 
be referred to as algebraic boson liquids. So far several 
algebraic liquid phases have been realized in toy models. 
In d = 3 stable algebraic spin liquids with gapless photon 
excitations have been realized in dimer model on a cubic 
lattice 2 , as well as spin models on a pyrochlore lattice 3 - 4 . 
In two-dimensional bosonic models with strong ring ex- 
change terms, an "excitionic bose liquid" phase has been 
proposed ir££. More recently another type of stable bo- 
son algebraic liquid phase with softened graviton exci- 
tations has been shown to be stable on an fee lattice^. 
Establishing the stability of such algebraic liquid phases 
is oftentimes quite subtle, since the gapless excitations 
are not due to the breaking of a continuous symmetry - 
they are not Goldstone modes^. 

The proliferation of topological defects in d > 1 often 
drives the instability of algebraic spin liquids. For exam- 
ple, as pointed out by Polyakov2ii£ ; the monopole (in- 
stanton) excitation in compact electrodynamics in d = 2 
will always proliferate and will gap out the photon exci- 
tation. On the other hand, if a bosonic model has a set 
of emergent symmetries at low energies which leads to a 
quasi-one dimensional behavior, as in the excitonic Bose 



liquid studied ir££, stability can be achieved by analogy 
with Luttinger liquids in d = 1. In the present paper 
we explore a new algebraic boson liquid phase in which 
the gaplessness is protected by a similar mechanism of 
an emergent quasi one-dimensional behavior. 

The model under consideration is motivated by cold 
bosonic atoms systems which can be trapped in square or 
cubic lattices formed by laser beams. The single particle 
ground state wavefunction within each well of the opti- 
cal lattice will have an approximate s— wave symmetry, 
with wave function <fio(r) ~ exp(— or 2 ). The first excited 
states are three fold degenerate p— wave states, with wave 
functions 4> x {f) ~ xexp(— or 2 ), 4> y {r) ~ ycxp(— ar 2 ) 
and 4>z{t) ~ z exp(— or 2 ). The wave function <f> x ex- 
tends further in the x than in the y and z directions, 
and so will preferentially hop to adjacent wells along the 
a;— axis. As discussed in**, although the p— wave states 
are not the ground states, the life time for these states 
can be much longer than the average tunnelling time be- 
tween neighboring optical lattice sites. If all of the atoms 
are pumped from the ground state to the p— wave states, 
they will establish a metastable equilibrium which can 
survive for considerable time. Here, following**, we ex- 
plore the properties of bosonic atoms trapped in such 
excited p— wave states. 

In previous studies, several phases have been identi- 
fied for this system, including a superfluid phase with 
1 dimensional I2 gauge symmetry 1 ^, the Mott insulator 
phase, a novel stripe phased 2 -. A supersolid phase has 
been predicted in even higher excited bands in cold atom 
systems in optical lattices 13 . In our work, the same sys- 
tem with strong interaction is revisited. The interaction 
induces a big gap to the excitations which change the 
total number of bosons on each site, and hence precludes 
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superfluidity. We establish that in addition to the Mott 
insulating state, a novel bond-algebraic liquid phase can 
be realized in both d — 2 and d = 3 systems. In these 
phases, bond-bond correlation falls off algebraically, but 
the correlation function is very anisotropic. We examine 
the instabilities of this gapless bond liquid phase towards 
both Mott Insulator and the Z2 phase. When the average 
filling per site is not a multiple of the spatial dimensions 
d, we find that, over a range of parameters, the bond 
liquid phase is stable. In d — 2, this liquid phase has a 
nice self-dual structure, which closely resembles the exci- 
ton bosc liquid phase studied earlier 5 . In the d = 3 case, 
we obtain a dual representation in terms of vortex loops 
which can hop on the dual lattice. 

This paper is organized as follows. In section II, 
the low energy effective Hamiltonian is derived, and the 
quasilocal symmetry of the low energy Hamiltonian is dis- 
cussed. In section III, we implement a spin- wave analysis 
which enables us to access and explore the bond algebraic 
liquid phase. In section IV, dual representations are de- 
rived in both d = 2 and d = 3. Section V is devoted 
to an analysis of the stability of the bond liquid phase. 
Instabilities towards both the Z2 phase and the Mott 
Insulator are also discussed. Employing the dual repre- 
sentation we identify the crystalline symmetry breaking 
pattern of the Mott insulator. In section VI, we consider 
modifying the model by doping with s— wave particles 
and softening the constraint on the particle number fluc- 
tuations within each well. Several additional phases are 
thereby obtained. Finally, section VII is devoted to a 
discussion of the prospects for experimental detection of 
the algebraic bond liquid phase. 



II. MODELS AND SYMMETRIES 

The Hamiltonian describing the cold atom system in- 
cludes the kinetic term, and the s— wave scattering: 



H= / d 3 x^(x)(- — V 2 - n + V T (x))il>(x) 
' 2m 



1 4na s h 2 

2 m 



d 3 x^(x)^(x)ip(x)ip(x), (1) 



where /u, is the chemical potential which determines the 
average filling of atoms n, Vt(x) is the optical trap po- 
tential, and a s is the s— wave scattering length. We now 
assume that within each well all of the atoms are main- 
tained in the p— wave excited states. One can then ex- 
pand the atom field operators in terms of the atomic 
p— wave states on each site 



(2) 



with a, b = x, y in d — 2 system, and a, b — x,y,z in 
d = 3 system. Here d- denotes a boson destruction op- 
erator of flavor a on site i, and i x ,iy, iz denote the x, y, z 
coordinates of the lattice site i. Inserting this expansion 



into equation ([T]), the kinetic term generates an on-site 
energy term for each flavor of atoms: 

(a)tj(o) 



<„ - l d 3 xMx)(~^- + VT(x))Mx)- (3) 

as well as a (dominant) nearest neighbor hopping term 
of the form, 



d 3 x4> a (x)( 



i a.b.c 

ft 2 V 2 



2m 



+ V T {x))(j> b {x-c). (4) 



The on-site energy for each flavor of atom needs not be 
the same, if the lattice symmetry is broken. For instance, 
if the laser beams in x and y directions have different 
wavelengths or amplitudes, p x particles should have a 
different energy from the p y particles. The hopping in 
(|U) is preferentially along one of the three axes, that is 
the integral in (|4|) is small unless a = b = c. The kinetic 
term is dominated by the very anisotropic hopping terms, 
for instance p x particles can only hop easily in the x 
direction. Since the other hopping terms will be small, 
we expect that over a large range of energies they will be 
unimportant. So we henceforth set these other hopping 
terms to zero. This is the central assumption of all that 
follows. Thus the hopping term is approximated as 



(«)t» 



h.c. 



(5) 



Due to the form of the p— wave orbitals the remaining 
single hopping strength, t' should be positive. 

In addition to the kinetic energy, the (s— wave scat- 
tering) interaction between the atoms will generate an 
on-site repulsion, as well as an intra well pair-tunnelling 
term between different flavors, 

H s = - 47ras? * / d 3 x^UxUUxU(xU(x) 
2 m J 

-EE M°M 6) + E 7(<4 a) M Q) M b ] 4 ] + h.c) 

(6) 

The U a b term encodes the repulsion between particles 
with the same flavor as well as between different fla- 
vors. The 7 term is the flavor mixing term, where, 
for example, two p x particles can be converted into two 
p y particles, and vice versa. The flavor mixing is due 
to the fact that, despite the vanishing overlap between 
the wave functions 4> x (t) and 4> y {r), the overlap integral 
/ d 3 x(4> x (r)*) 2 ((j)y(r)) 2 does not vanish. Angular mo- 
mentum conservation does not rule out this pair conver- 
sion term because, while the p x state is an eigenstate of 
L x , p y is not. Moreover, single particle flavor conversion 
is strictly forbidden being odd under parity symmetry, 

Px.4> x — ► -<b x , P x 6„ — > </>,,, 
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Py<f)y — > ~<Py, Py(f> x ~ > fa. (7) 

Because the integral J d 3 x(fa(r)*) 2 (4> y (r)) 2 is positive, 
7 will be positive. 

The interflavor and the intraflavor interactions will 
generally not be the same. For simplicity, we assume 
that the interaction matrix U a b together with the chemi- 
cal potential can be adequately approximated by a term 
of the form, 

a a 

^2fi a =n. (8) 

a 

Here n a is the average filling of each flavor. If all flavors 
are symmetric, n a = n/d. Also, let us not forget that 
since the wave functions of p a particles have compara- 
tively large overlap between two nearest neighbor sites in 
the same d axis (a = x, y, z), the s-wave scattering can 
generate a small off-site repulsion interaction between p a 

particles along d direction SH = a Sunf" n-'v^ , this 
anisotropic nearest neighbor repulsion is supposed to be 
much stronger than the nearest neighbor repulsion be- 
tween s-wave particles, which is usually negligible. 

Due to the flavor mixing 7 term in ((SJ) , the boson num- 
ber n x — n v is not conserved. However, as we will show 
below, in the phase of primary interest, the algebraic 
bond liquid phase, the flavor mixing term can be irrele- 
vant, therefore in this phase the boson numbers of each 
flavor of bosons can be considered as independent con- 
served quantities, and each flavor of boson has its own 
separate average filling. 

Let us assume U is the biggest energy scale in ([5]), 
and u is a small energy penalty for density fluctuation 
around average filling of each flavor, then the main effect 
of the term in is to keep the number of atoms in each 
well close to the average filling n, which we will always 
take as an integer. Moreover, we consider the strong 
correlation limit, where U is much larger than the kinetic 
energy term, t' . Since the average filling for each flavor 
is n a , when h a is not an integer, the system is effectively 
"fractionally" filled. Below, we will first focus on the case 
with h > 2, deferring a discussion of the behavior in the 
case with n — I, which is special because the intrawell 
pair conversion term alone cannot directly operate on the 
low energy subspace. 

The full Hubbard type Hamiltonian without SH then 
takes the form, 

H = if /, + Hu + H s = 

+ e - h 2 + E E 7(K (a)t ) 2 (rfl h) ) 2 + h.c.).(9) 

a i a^b 

In the large U = 00 limit, there will be a hard con- 
straint on the Hilbert space, with each site occupied by 



exactly n particles if n is integer. For large but finite U 
we will project the Hubbard Hamiltonian into this con- 
strained Hilbert space. A single particle hopping event 
will take one out of this low energy Hilbert space. At 
the second order, a "trivial" process which hops a p a 
particle at site i to site i + d, and then hops back to i 
can take place. This hopping process does not change 
the distribution of particles, but two p a particles at two 
nearest neighbor sites along a direction can benefit from 
this "trivial" process. Therefore this hopping generates 
a small effective attractive interaction between two par- 
ticles SH' = - J2i, a Svn^nj^, Sv ~ t' 2 /U. This off-site 
attractive interaction competes against SH, and gener- 
ally drives the system into orbitally ordered phase, if 5H 
is dominated by SH' , ferromagnetic distribution of fla- 
vors is favored; otherwise if SH' wins, p a particle at site 
i does not want to see the same flavor of particle stay at 
site i + d, in 2 dimensional case the system will be driven 
to the antiferromagnetic order of n x — n y . 

From now on let us focus on the case with Su w Sv, 
thus there is no obvious tendency to orbital order, the 
physics is controlled by higher order perturbations (our 
main result does not require Su exactly equals to Sv, we 
will argue that a small residual off-site density interaction 
does not destabilize the phase of main interest) . At forth 
order perturbation, a flavor-mixing ring exchange term 
of the form, 

Hring — 

-tYYid^dm^ J a) U {b) ^dfl & d^\ , + h.c.) 

Z_^ K 1 i+a i+a+b i+b i+b 1 l + a i+a+b ; 

(10) 

will survive the low energy projection. Here, t ~ 
(t') 4 /U 3 . The ring exchange in the XY plane is depicted 
in Fig. [TJ This term looks fairly complicated, but if we 
make the usual approximation, replacing the lattice bo- 
son operators with rotors d^ ~ e l6a , it will appear much 
simpler. The effective rotor Hamiltonian reads, 

H r = -t E E cos (9 a d b (6 ai - 9u)) + E u ( n ai ~ ^a) 2 

+ E7cos(2(0 Qi -0 6i ))(ll) 

a^b 

All the derivatives in (fTTj) should be understood as lattice 
derivative, i.e. da0(x) — 6{x + d) — 6{x). 

With n > 2, the Hamiltonian H r in equation (fTTj) is the 
starting point of our subsequent analysis. We first con- 
sider the symmetries of this Hamiltonian, initially tuning 
7 to be zero (the legitimacy of doing this will be discussed 
later). The symmetry in (jlip can then be summarized as 
follows. When d — 2, the Hamiltonian is invariant under 
the transformation, 

X -> Ox + h (y) + 9i (x) + a(x, y) 
V -» V + f 2 {x)+g2(y) + a{x,y), (12) 
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FIG. 1: A pictorial representation of the ring exchange term 
in (|1 If) . If a p x particle hops to the right, another p x particle 
must hop back to left, and similarly for the p y particles. 

and for d = 3 an analogous form, 

X ^>6 X + fi (y, z) + g x (as) + a(x,y,z) 
y -> 9 y + f%(x,z) +92{y) + a(x,y, z) 
0z^0z + h{x, y) + g 3 (z) + a(x, y, z). (13) 

Here fi,gi and a are arbitrary functions of their argu- 
ments. We can understand these symmetries as follows: 
the fi are functions of d — 1 coordinates reflecting the 
fact that each flavor can only hop in one dimension. The 
conservation laws associated with these U(l) symmetries 
correspond, for example, to the conserved total number 
of x— flavor particles at each fixed value of y and z. By 
contrast, the symmetries associated with the gi func- 
tions are additional emergent symmetries, present only 
after projecting into the low energy Hilbert space with a 
fixed number of particles on each site. Physically, in the 
projected Hilbert space each flavor of boson cannot hop 
freely: When one p x particle hops in the +x direction, 
another p x particle must hop in the opposite direction 
(see Fig 1). Finally, the a functions, depending on all d 
of the spatial coordinates, are (unphysical) gauge sym- 
metries corresponding to the hard projection imposed on 
the system, ^ a n a — n. 

If the 7 term is present, the continuous symmetries in 
(fl"2)) and (fl~3|) are broken down to Z2 symmetries. All the 
functions fi and gi can only take two values, or 71". 



III. ALGEBRAIC PHASES 

We first analyze the Hamiltonian in (fTTj) in the special 
case with 7 = 0. As we shall see, this enables us to access 
a novel algebraic bose liquid phase. We will then turn to 
the effects of a small but non-zero 7, and find that under 
certain conditions 7 is an irrelevant perturbation. Under 
a coarse graining in space and time it will scale to zero. 

To this end, we implement a "spin-wave" approxima- 
tion, expanding the cosine terms in (|11|) . The resulting 
action describes a Gaussian theory, with a Euclidian La- 



grangian of the form, 

L d=2 = ^{drOaf ^(d x d y (6 x - 6y))\ 



L d=3 = £ Y (d T 6 a ) 2 - £ - (d a d b (0 a - 6 b )) 2 . (14) 

a=x a=jtb 

Here we have rescaled space and time to set the velocity 
to unity. There is then just one remaining dimcnsionlcss 
parameter, K . Such a quadratic Lagrangian in which all 
terms involve derivatives of the fields, describes a scale 
invariant phase at long lengthscales. In a sense it can 
be viewed as a "fixed point" Lagrangian. The resulting 
phase is the desired algebraic boson liquid. 

The gapless modes of this Guassian theory can be read- 
ily calculated. In d = 2 there are two modes at each 
momenta with frequencies given by, 

<4 = 0, 

J( ~ sin(fc x /2) 2 sin(fc y /2) 2 . (15) 

The first nondispersing mode ujq is an unphysical gauge 
mode, corresponding to the function a(x, y) in (| L 2[) . The 
second mode vanishes along both coordinate axes in mo- 
mentum space, as required by the symmetries encoded 
in the functions fi and gi in (|12[) . This mode gives rise 
to a quasi one dimensional behavior, which is ultimately 
responsible for the stability of the novel bond algebraic 
liquid phase. If uj\ remains gapless in the presence of all 
allowed perturbations, this algebraic liquid phase is sta- 
ble. The double intraflavor conversion term and a set of 
"vertex" operators present in a dual representation (see 
below) can potentially destabilize the algebraic phase and 
gap out the mode uj\ . The stability of the algebraic liquid 
phase will be discussed later. 

In d = 3, in addition to the unphysical mode, uj® = 0, 
two other modes are obtained at each momenta. At small 
momentum k the dispersion relations take the form, 

lu\ ~ 2{k 2 x k 2 y + k 2 k 2 z + k 2 k 2 ), 

3 k 2 
2 _ 1.2 £.2 7,2 

W 2 ~ cy K x K y K 7. ,21.2 1 1,2 1,2 i 2,2 1.2 " 
z r^x^y ' ^x^z "r ^y^z 

(16) 

Notice that lu± vanishes along every coordinate axis in 
momentum space, whereas 0J2 vanishes on every coordi- 
nate plane. It is 0J2 that encodes the quasi-ld behavior, 
and as in d — 2 will be the ultimate reason for the sta- 
bility of the algebraic liquid phase. 

Due to the symmetries in (fT2"j) and (fT3"|) , many correla- 
tion functions vanish at this Gaussian fixed point. For in- 
stance, the correlation function between boson operators 
rf( x )t and <fW vanishes, as it breaks these symmetries. 
For the same reason, the correlation functions between 
two pair conversion operators vanish. But this does not 
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FIG. 2: A pictorial representation of the bond operator B x . 
If on one site there is one p x particle and on its adjacent 
site there is one p y particle, B x will exchange the p x and p y 
particles on the two sites. 



imply that this operator is irrelevant, because two adja- 
cent pair conversion operators can generate potentially 
relevant operators, which will be discussed later. 

The only correlators which do not vanish are correla- 
tion functions between bond operators. In the d = 2 case, 
the bond operators are defined as B x = cos(d x (9 x — 9 y )) 
and By = cos(d y (9 x — 9 y )). In the d — 3 case, the 
bond operators are defined in a similar way, for instance, 
B 



xy.x 



cos(d x (9 x — 9 y )). If on one site there is one p x 
particle and on its adjacent site there is one p y particle, 
B x will exchange the p x and p y particles on the two sites 
(Fig.©. 

The correlation functions between bond operators are 



d = 2 : {B x (0,0,0)B x ( T ,0,y)) 



1 



d = 3: (5^(0, 0,0,0)5 xy.x 

(r,0,j/,0)) 

1 



\ y 2 +T 2^'/K- 



(17) 



Here, the dimensionless number r]' depends on the form 
of the fixed point action on the lattice scale, that is not 
just the small momentum form . Notice that the correla- 
tion functions between two B x (B xy ^ x ) operators are only 
nonzero when they are at the same (r, y) plane. 



IV. DUAL REPRESENTATION 
A. d = 2 

In the previous section the Gaussian fixed point was 
obtained by expanding the cosine terms and allowing the 
fields to take on any real value. Here we address the 
legitimacy of this procedure. Specifically, we must take 
into account that the 9 fields are really phase angles (ie. 
compactified bosons) with the identification 9 = 9 + 2irZ. 
This implies that particular topological defects will be 
allowed. The effect of such defects can be most conve- 
niently addressed by passing to a dual representation, 
wherein one can identify "vertex" operators which insert 



such topological defects. The dual transformation of the 
theory should solve the constraint J2 a "-i^ = hi d = 2 
we define dual variables 6 and N as follows, 



d x d y ( 



ri, 

hi 



d x d y 4>, 
-d x d y </), 



(18) 



so that the constraint is automatically satisfied. Here, 
<f> and N are both defined on the dual site (plaquette 
centers). We can check the commutator and see that N 
and 4> are a pair of conjugate variables. As defined <fi takes 
on only discrete values, whereas N lives on the interval 
to 27r. But we will henceforth allow <fi (and N) to 
roam over all the reals, letting cj> — > ip and approximately 
imposing the discreteness by the addition of terms in the 
action. 

In terms of such a "coarse grained" field, ip, the Gaus- 
sian part of the dual action reads 



L d=2 = ^(dr^) 2 + J^(dxd y <p) 2 . 



(19) 



When n is integer, and both n x and n y are integers, the 
4>i should all be integers. Therefore in the dual Hamil- 
tonian there is a vertex operator cos(27riy9) on each site. 
However, since the dual Gaussian action (fH))) has a sim- 
ilar symmetry as in (|12[) . 

(20) 



¥> -> v + f(x) + g(y), 



the correlation function between two such cos(2nip) op- 
erators on different sites is zero. The vertex operators 
which show up in the low energy theory are dipole op- 
erators V x — cos(2ird x ip) and V y = cos(2ird y tp). At the 
Gaussian fixed point, the correlation function between 
two dipole operators is 



K(0,0,0)^(r,0,y)) 



1 



| r 2 +2/ 2|4K- 



(21) 



When h x and n y are not integers, <pi cannot all be 
integers. For simplicity, let us assume n is even and, 
n x — n/2+l/m and n y — n/2— 1/m, with m taking on all 
the integers. The arrangement of <f> can be taken as (j> = 
Z + xy/m. Now the potentially relevant vertex operator 
in the dual action is cos(2%d x cp+B x ). Here B x is a Berry's 
phase, B x — 2iry/m. When m > 1, this Berry's phase 
gives rise to an oscillating sign on the lattice, hence the 
vertex operator 008(2^9x^5 + B x ) does not appear in the 
low energy coarse grained theory. However, vertex oper- 
ator Vq PtX = cos(2qnd px ip) = cos(2qirip(p, 0) — 2qirip(0, 0)) 
(with integer p,q), as well as "parallel jump" (PL)term 
V p i = cos(27r9 x (ya(0, 0) — 2Trd x <p(l, 0)) do not contain any 
oscillating sign on the lattice as long as qp is a multiple 
of m, and hence should appear in the low energy theory. 

The leading (most relevant) operators are V\ m .x arid 
V m i. x , as well as V p i operators. The correlation functions 
between V qp , x and V p i are 



(V qp , x (0,0,0)V qPtX (T,0,y)) 



1 



<Vpi(0,0,0)Vpi(T,0,i/)) 



1 



! r 2 _|_ y2|32if/3 ' 

Here 8(j>) is an integral which depends on p, 

(2sm(pk x /2)) 2 



1 



S (p) = o / dfc * 



sin(fe x /2) 



(22) 



(23) 



When p is large, <5(p) roughly scales as S(p) ~ lnp, there- 
fore, in the case of large p, V p i is the most relevant oper- 
ator in the dual formalism. 



FIG. 3: The dual theory defined in (|24p describes vortex loops 
which can only hop across a cube, as shown above. 



B. d = 3 

In the d = 3 case, the compactification of the phase 
angles can also be dealt with in a dual representation. In 
the case with d — 3, we only focus on the situation with 
n x — n y = n z — n/3, as in this case, so long as n ^ 3Z, 
it is enough to guarantee a stable liquid phase. Let us 
define dual variables and 4>i a $ 

N z = d x d y (e x -e y ), 
N y = d z d x (e z -8 X ), 
N x = d y d z (9 y -e z ), 

d z d x (j>y - d z d v (j) x = n/3 - n z , 
d x d v cj) z - d x d z (j)y = n/3 - n x , 
d y dz4>x - dyd x (f) z = n/3 - n y , 

[cl) i ,N j }=i6 ij . (24) 

Here N x , N y and N z can be interpreted as elements 
of "vortex" lines in the 9 a configurations, along the x, y 
and z directions, respectively. Moreover, cos(d z d x (j) y — 
dzd y 4> x ) is a term which hops a unit size vortex loop in 
the xy— plane along the z direction (Fig. [3]). Therefore 
the dual picture describes a system of vortex-like loops 
which can hop in one dimension. This is analogous to 
the roton hopping model studied in d = 2^. 

Various phases can be understood in terms of vortex 
loops in boson systems. The normal three dimensional 
boson theory is dual to a theory of vortex loops minimally 
coupled to a rank-2 antisymmetric gauge field. The su- 
pcrfluid phase is a phase in which the vortex loops tend 
to shrink, and the gapless rank-2 antisymmetric gauge 
field corresponds to the Goldstone mode in the original 
theory. In the 3 dimensional crystal (or Mott insula- 
tor) phase of the bosons, the vortex loops proliferate and 
"condense" , gapping out the rank-2 antisymmetric gauge 
field through the Higgs mechanism. If there are two fla- 
vors of vortex loops, a gapless boson liquid phase can be 
realized, where the fluctuation of vortex loops become the 
photon excitations in the photon boson liquid 15 . Here we 
show that when the vortex loops hop in a perpendicular 
direction only, we can recover a dual representation of 
the algebraic boson liquid of the original p— wave boson 
system. 



The Gaussian part of the dual vortex loop hopping 
action reads: 



Lloop — 

(25) 

Again, ip a is a coarse grained field of <j) a . The symmetries 
of this dual Gaussian action are 

fx -> (fx + h(x,y) + f 2 (x,z) + d x a(x,y,z), 
(fi y -> tp y +gi {x,y) + g 2 (y, z) + d y a(x,y, z), 
¥>z -> 9z +hi(x,z) + h 2 (y,z) + d z a(x,y,z). (26) 

Due to its definition in (|2"1|) . <\> a only takes on dis- 
crete values. Thus, as in the d = 2 case, vertex oper- 
ators are allowed in the dual theory. Most of these ver- 
tex operators are irrelevant at the Gaussian fixed point 
since the correlation functions between them break the 
symmetries listed in (|26p. One important vertex op- 
erator with nonzero correlations is a loop growing op- 
erator Gab = cos{2ir(d a (pb — d ip a ) + £>), where B is 
again a Berry's phase. This loop growing term can en- 
large, shrink, annihilate and create a vortex loop (Fig. 
3]). A Berry's phase term will generally be present in 
the dual vertex operators, except when n = 3Z where 
it vanishes. Since the filling of the boson is n/3 for 
each flavor, when n is not an multiple of 3, the vertex 
operators which survive in the low energy field theory 
should be G 3 i ia h = cos^^a^b - d b ip a )) and Gi 3i0 b = 
]li=o cos(27r(<9 a </? b - d tp a )( a ^ +l )). Gi 3iQb is an operator 
that creates a "large" vortex loop which covers three unit 
squares (Fig. [6]) . Another potentially relevant vertex op- 
erator is the parallel jumping term, denoted as V p i. This 
term hops a unit vortex loop parallel to its unit square 
(Fig [5]) . Both the loop growing term and the V p i term 
can only correlate in the direction perpendicular to the 
loop. The correlation function for these two operators 
are, 

(G 3 i,^(0,0,0,0)G 3 i,^(r,0,0,^)} 
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FIG. 4: The loop growing vertex operator can enlarge a vortex 
loop. The arrow is the direction of the vortex. 




FIG. 5: The parallel jump vertex operator can jump a unit 
vortex loop to its adjacent unit square. The arrow shows the 
direction of the vortex. 



\t 2 +Z 2 \ k ^' 



(Gi 3 ,^(0,0,0,0)Gi3,^(t,0,0,z)) 

1 



T + Z 



2\K m 



CMo,o,o,o)-M^o,o,z)> 



\t 2 + z 2 \ k ^' 



(27) 



FIG. 6: The "large" loop creation operator Gra, a b, which cre- 
ates a "large" vortex loop covering three unit squares. 



Again, rji are coefficients which depend on the lattice 
scale physics. The smallest of the three corresponds to 
the most relevant operator in the dual formalism. The 
numerical result gives that 773 = 1 1. 3 is the smallest coef- 
ficient, therefore V p i is the most relevant operator in the 
dual formalism. 



V. INSTABILITIES AND RESULTANT PHASES 
A. Instability 

The algebraic phase has two dominant instabilites, 
which are expected when the dimcnsionlcss coupling con- 
stant K is very large or very small. If the pair tunnelling 
term in (jTTJ) is relevant, it breaks the symmetries listed 
in (fT3")l and (fT2"|) down to Z2 symmetries, which then gaps 
out the gapless excitations of the algebraic phases. The 
correlation functions between two pair-tunnelling terms 
7 cos(2(# a — 9 b )) vanishes within the algebraic bond liquid 
phase, but two adjacent pair-tunnelling terms can gener- 
ate pair-bond operators such as B 2x — cos(2d x (8 x — 6 y )) 
and B2ab,a = cos(2d a (0a, — 9b))- The correlation functions 
between pair-bond operators take the form, 



d = 2: (B 2x (0,0,0)B 2x (T,0,y)) 



1 



| r 2 +y 2|4/(7r 2 iTj> 

d = 3 : (B 2xyiX (0, 0, 0, 0)B 2xv , x (t, 0, y, 0)} 

1 

~ \ij 2 + t 2 \^'/k- 

(28) 

When the pair-bond terms become relevant and grow in 
magnitude under coarse graining, the 7 term itself will 
become important too, and the system is expected to 
enter a Z2 phase, which we discuss in the next subsection. 
The pair-bond terms are relevant when K > 2/ir 2 for 
d = 2, and K > 2r( for d = 3, so the Z 2 phase will occupy 
the large K region of the phase diagram as shown in Fig. 
6. 

The vertex operators in the dual representation can 
drive another instability when K is small. Once rel- 
evant, these vertex operators destabilize the algebraic 
phase, giving a Mott insulator which will break lattice 
symmetries due to the Berry's phases. The most rel- 
evant vertex operators in the dual representation de- 
pend on the dimension and the average filling. When 
d = 2 and n a are integers, the most relevant terms are 
V x = cos{2ird x ip) and V y = cos(27n9 y </?); when d = 2 and 
to is large enough, the most relevant terms are paral- 
lel jump terms V p i. If the vertex operators are relevant, 
they will obviously break the symmetry in (|2H)) . since if 
we expand these vertex operators, they would read 



V X - {d X ip) 2 ,Vy - (dyifi) 2 
V pl ~(d 2 M 2 + (d 2 y<P) 2 



(29) 
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Mott Insulator Bond Algebraic Phase Z 2 phase 



Featureless MI Z2 phase 

+ ► 

K c K 

FIG. 7: The phase diagrams for the p-wave particles in both 
d — 2 and d = 3 case. The top diagram is for the case with 
fix — n y = 2/m, and m > 1 in d = 2; or n 7^ 3Z in d = 3. In 
both cases nonzero Berry's phases enter the vertex operators, 
determining the crystalline order in the Mott insulator. The 
bottom diagram is for filling n = dZ, in which the Berry's 
phase vanishes and the Mott Insulator is featureless. Since 
the vertex operators are more relevant when n — dZ, the 
algebraic bond liquid does not exist as a stable phase. 

These terms do not have the one dimensional symmetry 
any more. 

When d = 3 and n — 3Z, the most relevant terms 
are loop growing terms G a b = cos(27T {d a tpb ~ dbtp a )), 
and when d = 3 and n ^ 3Z, the most relevant terms 
are again the V p i terms. Based on the correlation func- 
tions in (|22]) and (f27|) . the V p i terms are irrelevant when 
K > 3/16 for d = 2 and K > 2/773 for d = 3. The scaling 
dimensions of V\ m and G13 operators are higher than the 
V p i terms as long as m is big enough. A numerical eval- 
uation of the dimensionless numbers, which involves an 
integration over the full Briullioun zone, gives 773 = 11.3 
and rj' = 0.0950. Therefore there are finite parameter 
regions for the existence of a stable algebraic bond liquid 
phase in both d — 2 and d = 3 (Fig. [7|) : 

d = 2 : 3/16 < K < 2/vr 2 , 

d = 3 : 2/7/3 < K < 2t/. (30) 

However, when h — dl where the Berry's phase terms 
vanish, lower order vertex operators are allowed which 
are more relevant. In this case the bond algebraic liquid 
will not exist as a stable phase (Fig. [7]). One notable 
aspect of this theory is that the calculation of the scaling 
dimensions depends on the behavior of the model on the 
lattice scale. 



B. Z2 phase 

When the coupling K is sufficiently large the B 2x and 
B 2a b,a terms will be relevant, and the pair tunnelling 
term in (|11[) will also grow large and will drive the system 
into the Z2 phase. What is the nature of this phase? First 
we notice that, even if the pair conversion term is rele- 
vant, the symmetries in (fT2")l and (fT3"|) are still respected 



as long as /j and gi take only two values, or n. In the 
following, we will first focus on the d = 2 case, which is 
technically much simpler. One can directly generalize to 
the d = 3 system. 

Since 7 should be positive, when it becomes large it 
will lock 9 X — 9 y into two values, ±7r/2. Because n x + n y 
is a static variable in the projected Hilbert space, we can 
just focus on the variables n x — n y and 9 X — 9 y . Notice 
that if n x + n y is even (odd), n x — n y has to be even (odd) 
too. 

Let us assume n x + n y is even, for simplicity. Two Z2 
variables can be introduced to describe the system: 

a z = i exp(i(9 x — 9 y )), a x = exp(iir(n x — n y )/2). (31) 

Using the formula 

e A £ B = e [A,B] e B e A (32) 

we can see a x and a z anticommute with one another, so 
they have the same algebra as Pauli matrices. In terms 
of these "spin" variables, the original model (jTTJ) can be 
effectively written as, 

H z2 = £ -hat - E K<t°2<4°l- (33) 

i □ 

This model was studied previously in reference [la] , where 
it was introduced to describe a phase transition in an 
array of p + ip superconducting grains. For K > h, it 
was shown that the system is in a bond ordered phase, 
with the expectation value of a*a* +x being independent 
of Hi. This phase has a large degeneracy, D, which can 
be estimated by counting the number of Z2 symmetries 
of (HU). The result is 

D = 2 L*+Ly (34) 

The physical transformations within the degenerate 
ground state manifold only transform 9 X — 9 y . One of 
the classical ground states when h = is shown in Fig. 
[5|i. When K < h, there is no bond order, and the ground 
state is not degenerate. The transition between these two 
phases happens exactly at K = h. 

In d = 3 with 7 > 0, the three pair tunnelling operators 
are 

H v = 7(008(2(0., - By) + cos(2(0 y - 9 Z ) + cos(2(0 2 - 9 X )) 

(35) 

These three terms cannot be independently minimized, 
i.e. they are frustrated. Minimizing the sum of all three 
terms gives the following solutions, 

9 X - 9 y = 7rfc 3 /3, 9 y - 9 Z = irki/3, 9 Z - 9 X = irk 2 /3 1 

k 1 +k 2 + k ? , = 0,(36) 

where the ki are not divisible by 3. The symmetry relat- 
ing the manifold of solutions is the same as in (|13p except 
that fi + gi only takes two values, or tt. 
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FIG. 8: One of the classical ground states of model Q33J) when 
h — 0. The ground states satisfy the condition that each 
unit square on the lattice has either (2+, 2—), or (4+, 0—) or 
(0— ,4+) configurations. 



As in d = 2, this Z2 phase exhibits bond order, with 
spatially varying B a b, a and -B a fc,b- For instance, the ex- 
pectation value of B X y iX — cos(d x (9 x — 9 y ))(x,y, z) is 
independent of y. The expectation value can be ±1 or 
±1/2. Once again, this Z2 phase has a large degener- 
acy. In the Z2 version of the symmetry (|13p . the func- 
tions fi and gi correspond to physically different states 
in the ground state manifold, but this is not the case 
for the function a. The total degeneracy D is obtained 
by counting all the physical symmetry transformations, 
which is 



(37) 



When h — 1 the pair tunnelling term takes the system 
out of the projected Hilbert space. But acting in concert 
with hopping processes it will generate terms such as 



H e = jcos((9 x - 6 y )i - 

+(e x -9 y ) 3 + (6 



'V)2 



(38) 



Here 1,2,3 and 4 label the four corners of one unit square. 
This process is shown in Fig. [9] The correlation function 
between two H e operators is zero in the algebraic phase, 
but H e terms can generate terms which have nonzero 
correlation functions along certain directions. If we check 
carefully, we will find that the leading term generated is 
still irrelevant enough for a stable algebraic liquid phase. 
Thus, with h — 1 the algebraic boson liquid can exist as 
a stable phase. 

What kind phase will H e drive the system into when 
it is relevant? Notice that equation (f38|) has the same Z2 
symmetry as (|33[) . with the particle number of n x — n y 
conserved mod 2 along every row and column. Thus the 
order should again be bond order. Indeed, if we take 
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FIG. 9: The process in (|38[1 which involves four intermediate 
states. Step (a), p x particle at site 3 hops to site 4 (interme- 
diate step); step (&), two p x particles at site 4 can be trans- 
formed to two p y particles due to the pair tunnelling term; 
step (c), p y particle at site 2 hops to site 3; step (d), p x par- 
ticle at site 1 hops to site 2; step (e), one of the p y particles 
at site 4 hops to site 1). 



and the ring exchange term ([IT]) together, the Hamilto- 
nian involving 6 = 9 X — 8 y reads, 



H 



-t cos(6*i 
'cos(— 9\ 



9 4 ) +7Cos(6»i 
9 4 ) +7Cos(6»i 
+7Cos(6>i 



^2 1 (73 ^4 / 

h + O3 + 4 ) 
h - O3 + 04) 



(39) 



One of the minima of this Hamiltonian is 0j = n/2. All 
of the other minima can be obtained by transforming 
9i —> 9i + ir on each row and column. The symmetry and 
the ground state are thus the same as (|3"3"1) . The resulting 
phase is precisely the Z2 phase found before. 



C. Mott Insulator phase 

For small enough coupling K, when the vertex oper- 
ators in the dual formalism are relevant, the system is 
driven into the Mott insulator phase. In the integer fill- 
ing case with n = cK, the Mott insulators are featureless 
and there is a unique ground state with no broken sym- 
metries. With fractional filling, flavor-density ordered 
patterns which breaks space symmetries are expected. 
Again, let us first focus on d = 2 case, and then general- 
ize the results to d = 3. 

In the 2 dimensional Mott Insulator phase, because the 
average filling of each flavor of bosons could be fractional, 
crystalline order is expected. However, the crystalline 
order is only in the n x — n y channel, since n x + n y is 
fixed to be n. The dual action in the case with n = 2Z 
and n x — n y — 2/m reads, 



Lioop = J^idrV) 2 + j^(dxd y ip) 2 



cos(2mird x ip) — 7 cos(2rmtd y (p) 



(40) 



The ellipses include other vertex operators, for instance 
the V p i terms, however, these terms do not determine the 
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FIG. 10: The stripe order of n x 
odd integer and n x = n y . The 
with 71^ 



— n y in the case of n is 
full circles represent sites 
> 0, and the empty circles represent sites with 



n„ < 0. 



ordered pattern. The double vertex operators induce de- 
generate ground state once they are relevant. The most 
natural phases have stripe order in the n x — n y density 
or a plaquette order. The order parameters can be rep- 
resented in terms of low energy variables: 



x + n y )e 27rlx/m ~ sinpTrfl^), 



cos{d x d y {6 x - 6 y ))a 2 ™l m ~ cos(27r9^), 
cos(d x d y {9 x - 6 y ))e 2 ™ y l m ~ cos(27r9^). 



(41) 



This representation can be obtained from either sym- 
metry arguments or from renormalization group flow of 
the lattice theory. If the vertex operators in (|4"0")) are rel- 
evant, the stripe order parameters in (|4ip take nonzero 
values. If 7 in (|40|) is positive, the plaquette stripe order 
is favored; if 7 is negative, the n x — n y density stripe or- 
der is favored. If n is odd integer and n x = n y , the dual 
action and the low energy representation of order pa- 
rameter are similar to the equations (|40|) and (|4"Tj) with 
m = 2. The n x — n y density stripe order with odd inte- 
ger n is depicted in Fig. 1101 Both density stripe order 
and plaquette stripe order have been found numerically 
in a spin-1/2 model with a similar 1 dimensional gauge 
symmetr y i 19 . Notice that the crystalline pattern pre- 
dicted from the field theory is only valid when the cor- 
relation length is much longer than the lattice constant, 
i.e. one is close to the phase transition between the Mott 
and algebraic liquid phases. 

In d = 3, in the case of n 7^ 3Z and n x — n y = h z — 
n/3, the dual action including the vertex operators reads, 



L 



loop 



J2j^( d rVa) 2 + Yj ^( d a{db(Pc ~ d c ip b )) 



- y^gcos(67r(d a y? fc 



dbWa)) 



(42) 



a^b.b^c,c^a 



The triple vertex operators gives rise to degenerate 
ground state when they are relevant. Again, because 
n x + n y + n z is fixed at n, the density modulation can 
only happen in the n a — n& variables. The most natural 
state is the stripe order, and the order parameter can be 
represented in terms of low energy variables as, 



n x — n y ~ exp(27ria;/3) sin(27r(<9j,<y9 z — d z tp y )) 

- exp(27ra//3) sia.(2ir(d z (p x - d x tp z )), 

n y — n z ~ exp(27riy/3) sm(2ir(d z ip x - d x ip z )) 

- exp(2niz/3) sm(2n(d x (p y - d y ip x )), 

n z — n x ~ exp(27riz/3) s\n(2i:{d x ip y — d y ip x )) 

- exp(2nix/3) sm(2ir(d y tp z — 9 z (p y )), (43) 

The plaquette stripe order parameters can be obtained 
in a similar way. These order parameters take nonzero 
expectation value once the vertex operators in ()42|) are 
relevant. Because the order parameters in (|4"3"|) are not 
invariant under the space symmetry group, the ground 
states will have spontaneously broken symmetries. 

A small off-site density-density interaction between 
bosons can in principle exist, e.g. the SH and SH' terms 
discussed in the second section. In the bond liquid phase, 
a term of this form renormalizes the field theory param- 
eter K, but it does not introduce any new vertex opera- 
tors. Thus a perturbation of off-site density interaction 
does not destabilize the liquid phase. 



D. Other instabilities 

Throughout the paper, we have assumed that the hop- 
ping amplitudes, t a ^,c in (U), are only non-zero when 
a = b = c. In a real system, t a j,^ c will also be non-zero 
when this condition is not satisfied, but are expected to 
be very small. Once non-zero, these hopping processes 
will lead to new instabilities of the bond algebraic liq- 
uid phase. For instance, at the order of t XjXtX t yt y tX /U, 
bond operator B xyx — cos(d x (6 x — y )) is generated. 
Once present this term will be relevant in the bond al- 
gebraic liquid phase. The system will be driven into a 
phase where e l ( e *>~ Sy ^ = ±i, the ground state is two fold 
degenerate. But so long as these other hopping ampli- 
tudes are very small, this instability will only occur at 
very low temperatures, and at higher temperatures the 
system will behave as if it were in the algebraic liquid 
phase. For instance, at small temperature the system 
can be in the quantum critical region controlled by criti- 
cal point at which these perpendicular hopping terms are 
tuned to zero, this quantum critical point is described by 
the Gaussian actions (CO 
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VI. OTHER PHASES 
A. Doping with s wave particles 

The results in the previous sections have been obtained 
under the assumption that all the particles are in p— wave 
states. But what behavior is expected if some of the par- 
ticles are in the s— wave ground state? Let us take d = 2 
as an example. Let us suppose that the total lattice fill- 
ing is still an integer h. For strong repulsive interactions 
which preclude superfluidity, we again project into the 
constrained Hilbert space with ^ a n a + n s = n on ev- 
ery site of the lattice. Here n s is the density of s— wave 
particles, and we will denote the conjugate phase field as 

Since the s— wave particles can readily hop in all di- 
rections, the low energy effective hopping within the pro- 
jected Hilbert space will allow a term of the form (assum- 
ing d = 2 for simplicity): 

H h = ~h cos(d x s - d x x ) - t x cos(d y 9 s - d y 9 y ). (44) 

The presence of the s— wave particles can thus strongly 
modify the dynamics of the p— wave particles. Moreover, 
in addition to the allowed pair tunneling between the p x 
andpj, flavor particles, it will also be possible to exchange 
p— wave and s— wave particles. Although a single particle 
tunnelling is still forbidden by parity symmetry, pair tun- 
nelling between s— wave and p— wave particles will be al- 
lowed. Upon including these processes, the full Hubbard- 
like Hamiltonian in the rotor representation will take the 
form, 

H r = -tcos(d x 9 s - d x 9 x ) ~ tcos(d y 9 s - d y 9 y ) 
y 

} J u p (n a - n p /2) 2 + u s (n s - n s ) 2 

a—x 

V 

~^g a cos(2(0 s - 9 a )) - g xy cos(2(0 B - 9 y )) (45) 

a—x 

Here n s and n p are the total fillings for the s— wave and 
p— wave particles, with n s + h p = h. We are implicitly 
assuming that n > 1. 

Consider first the situation in which the s— wave par- 
ticles are at integer filling and have a charge gap, which 
will be the case for large enough u s . The above hopping 
term will then take the system out of the low energy 
sector, but can generate an effective ring exchange term, 

H r i n g — —t cos(d x d y (9 x — 9 y )) -\ . (46) 

This term is the same as the ring exchange term in (|lip . 
Therefore when the s— wave particles are gapped, the sys- 
tem behaves similarly to the undoped expected 
on physical grounds. 

To address the more interesting situation in which 
n s ,n p are non- integer, it will be helpful to first consider 
the symmetries of the above Hamiltonian in the absence 



of the pair tunneling terms, g a — g xy — 0. In this case 
the model is invariant under, 

9 X -> 9 X + fi(y)+g%(x) + a(x,y), 
9 y -^9 y + f 2 (x) + g 2 (y) + a(x, y), 
9 s ^9 s + gi (x)+g 2 {y)+a(x,y). (47) 

We now search for a stable algebraic phase by first set- 
ting g a = g xy — 0, and then expanding the cosine terms 
in H r to obtain the following quadratic Lagrangian, 

L d = £ |(d T ^) 2 + ]Tf (d Q (0 s -0 Q )) 2 . (48) 

fj,=x,y,s a—x 

Upon diagonalization, one obtains 3 modes, which at 
small k are given by, 

^o 2 - 0, 

k 2 k 2 

w i ~ ^x + Kj ~ \J k i ~ k%ky + ky- ~ ^ 2 , 

lo 2 ~k 2 x + k 2 y + Jki -klkl + k y ~ k 2 . (49) 

Again luq corresponds to the unphysical gauge mode re- 
sulting from the hard constraint ^2 a n a + n s — n. The 
u>2 mode corresponds to the gapless Goldstone mode re- 
sulting from the "condensation" of 9 S . The uj\ mode 
vanishes along the coordinate axis in momentum space, 
and it gives rise to quasi-ld physics. 

The relevance of the pair tunnelling terms can be es- 
timated by considering their scaling dimension in the 
Gaussian phase action ([48]) . Based on the above sym- 
metries of the Gaussian theory, it is concluded that the 
operators cos(2(# s — 9 a )) can only correlate in one di- 
mension. Moreover, a calculation shows that it has a 
power law correlation function, due to the contribution 
from uj\ in (|49[) . The parallel jump term generated by 
cos(2(# Q — 9 )) is not as relevant as cos(2(0 s — 6 a )). When 
cos(2(# s — 6 a )) is irrelevant, the system is in another type 
of algebraic liquid phase different from the bond liquid 
phase, which is described by the Gaussian action (I48|) . 
The operators which have power law correlations are the 
tunneling operators cos(2(0 s — 9 a )) instead of bond op- 
erators. 



B. Id liquid phase and superfluid phase 

So far everything we have considered is under the con- 
straint ^ a n a — h on every site of the lattice. What 
will happen if we soften this constraint? Softening the 
constraint will lower the gap for the charged excitations, 
and eventually superfluid order can develop which breaks 
the global U(l) symmetry of the phase J2 a @a- This su- 
perfluid phase has been discussed in reference^. But as 
we now show it is also possible to have a stable quasi 
one-dimensional liquid phase. 



12 



Let us take the d = 2 system as an illustrative example. 
The original Hamiltonian before projection is © . We 
notice that, with the pair tunnelling term, the symmetry 
of the system is the Z2 version of the Id symmetry, and 
there is moreover a global 2 dimensional U(l) symmetry 
Ox — * x + 9q, 9y — * 0y + do, for a spatially independent 9q. 
In the case with small interaction U, upon passing to the 
rotor representation, and expanding the cosine hopping 
terms, the action associated with the Hamiltonian in ([9]) 
becomes, 

a— x,y 

-7cos(20 x -20„). (50) 

When 7 = this Lagrangian describes a type of one- 
dimensional liquid phase, which can be visualized as a 
set of decoupled 1 dimensional Luttinger liquids on every 
row and column of the 2 dimensional square lattice. To 
examine the stability of these Luttinger liquids one must 
check to see when the vortex tunneling events in the dual 
representation are irrelevant^. If we assume n is odd, 
the average filling for each flavor is half-integer, and the 
single vortex operator in the dual theory has a nonzero 
Berry's phase. Double strength vortices are irrelevant 
provided 20 K < 2ir. For K > 2tt the proliferation of 
these vortices will destabilize the quasi Id liquid phase 
and drive the system into a Mott insulator with broken 
translational symmetry. 

Stability of this quasi-ld liquid phase also requires the 
irrelevance of 7. First we note that the two-point cor- 
relation function of the operator cos(20a; — 26 y ), when 
evaluated in the quadratic theory, vanishes unless the 
operators are on the same site. We can thus focus exclu- 
sively on the time dependence of this correlator at (any) 
one given point on the 2d lattice, say x = y = 0. This site 
lies on only two of the Luttinger liquids, one running in 
the x— direction and the other in the y— direction. This 
problem can then be viewed as a zero dimensional point 
contact between two Luttinger liquids, which has been 
studied in detail. Following this work we find that 7 is 
irrelevant when K > ir. We thus conclude that when n 
is odd, there will be a range of parameters, n < K < 2n, 
within which the quasi-ld liquid exists as a stable phase. 

What if 7 is relevant? In this case we should mini- 
mize the cosine potential, which will lock the two fields 
+ run/2. r)i = ±1, and can be represented as 
r/i = (— l)/(^)+f(' i «). f and g are both integer valued 
functions. In the simplest case wherein / and g are both 
constants, then the quadratic part of the Lagrangian 
when re-expressed in terms of a new field, = 6 x +6 y , will 
be simply L ~ Ylu=x y t(^m^) 2 - This is the Lagrangian 
of a conventional 2d superfluid, found in previous worbii. 
We thus conclude that the 7 term tends to drive the sys- 
tem from the quasi-ld liquid phase into a 2 dimensional 
superfluid phase. 



VII. EXPERIMENTAL IMPLICATIONS 

Here we discuss how one might try to detect experi- 
mentally the presence of the algebraic bond liquid phase. 
We will consider three different types of correlation func- 
tions which can in principle be measured in cold atom 
experiments. First, we will consider the particle density- 
density correlator, which can in principle be extracted 
via light scattering experiments. Secondly, we will con- 
sider the momentum distribution function, = (ipt'fpk)- 
In cold atom experiments this quantity can be extracted 
by releasing the atoms from the trap, and at some later 
time measuring their spatial distribution function in real 
space, n(R), with R the radial coordinate measured from 
the center of the trap. At time t after releasing the atoms, 
these two quantities are related as, nk= m R/t = n(R). Fi- 
nally, we consider the momentum density - momentum 
density two point function, (riknk'}, which can be mea- 
sured by spatial noise correlation o 21 ' 22 . 

A. Particle density-density correlation function 

The particle density-density correlation function 
(p(O)p(r)) can be used to distinguish between various 
phases. For example, the crystalline Mott insulator 
phases manifest themselves as Bragg peaks, whereas 
quasi-ld liquid phase results in a power law behavior. 
The Bragg scattering techniques have been applied to 
atoms trapped in an optical lattice 2 ^, and the periodic- 
ity of optical lattice was extracted. The particle density- 
density correlation function depends on the form factor 
of p— wave particles, and the flavor density-density corre- 
lation function (for instance (n x (i)n y (j))). In the phase 
of primary interest - the algebraic liquid phase - the cor- 
relator between total flavor density operators provides no 
useful information since the total particle density on each 
site is fixed at n. Thus, the correlation functions between 
the total density operator n x + n y will be quite trivial. 
However, since p x and p y particles have different form 
factors, the particle density-density correlation function 
(p(0), p(r)) contains nontrivial information about this al- 
gebraic phase. 

We first briefly consider here the behavior of the corre- 
lator of the flavor-density operator, which is non-trivial 
in the algrebraic liquid phase. From the d — 2 duality 
we have, 

n x -n x = d x d y 4>, 
n y - n y = -d x d y (j>. (51) 

In terms of the coarse grained field ip, the representation 
of n x and n y are 

n x ~ n x + d x d y p 
+c(e 2vix / m sm(2TTd y p) + e 2my / m sm(2nd x p)), 

n y ~ fiy- d x d y p 
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z (e 2nx l m sin(27r<9^) + e 27riv/m sm(2iTd x cp)). 



(52) 



The correlation functions can then be readily computed 
giving, 

((n x (0, 0, 0) - n x )(n x (x, 0, r) - n x )) 

2-Kix/m 



In the case of Mott Insulator phase with a crystalline 
pattern of flavor difference n x — n y , the result is 



(Pk,P- 



S(k - G[){kl - k 2 y ) 2 + c 2 S(k G t ). 



(56) 



Here G\ are the crystalline order wave vectors. Therefore 
we conclude that using light scattering experiment, one 
can detect the bond algebraic phase. 



((na. (0,0,0) - n x )(n x (0,y,r) - n x )) 

2iviy/m i „ _ 



r 2\4K 



(M0, 0, 0) - n/2)(n x (x, y, 0) - ft/2)) 



X 2y2 



(53) 



When the vector separating the two n x operators is 
parallel to one of axes, x or y, the correlation function 
has two contributions. The first comes from the sine 
term, the second from the derivative term in (|52j) . 

When the vector separating the two n x operators is 
not along an axis, only the derivative term in (|52|) con- 
tributes. Notice that n x not only correlates in the x di- 
rection, but also in the y direction. This is due to the fact 
that, the ring exchange term in (|11[) has no asymmetry be- 
tween x and y. However, in the correlation function (|53p 
the two dimensional rotational symmetry is not restored 
even when the scale is much larger than the lattice con- 
stant, which is quite different from other algebraic liquid 
phases. 

The particle density correlation functions are supposed 
to be measurable with light scattering experiment. The 
result is usually written in momentum space (pk,P-k), 
which contains information of both the form factor of 
p— level state wave functions, and the particle density 
correlation functions. In the case of the bond- algebraic 
phase, the result at small momentum is 

(Pk,P-k) 

-EE («*.«#) x / d 2 ^V^-%»(^)l% fc (y)l 2 

i,j a, 6=33 

~ F(k x , k y ) x {k 2 x - k 2 ) 2 + J2 c i S (k - Gi). 



B. Momentum distribution function 

The momentum distribution function has been used to 
detect both the Bose-Einstein condensate, being charac- 
terized by a sharp peak at k — 0, and the Mott insulator 
where there is no sharp peal*22. We can calculate the 
momentum distribution function of various phases as fol- 
lows. Focussing for simplicity on 2 dimensional systems, 
the momentum distribution for the superfluid phase of 
s— wave particles is 

n{k)/N=l/N { d 2 r 1 d 2 r2(ipHri)^(r2))e l '* iPl ~ p2) 



~ 1/N J d 2 r 1 d 2 r 2 exp(iA: • (fi - f 2 )) 

xE(^(^)4 s)t ^(^)4 s) ) 

l/7VEe^^-^^ s (fe)| 2 (4 S)t 4 S) > 

~Y,\Mk)\ 2 6(k-G t ) (57) 



i,3 



(54) 



Here Gi are the vectors in the reciprocal lattice of the 
optical trap potential, and F(k x ,k y ) is a function of k x 
and k y , 

y 

F(k x ,k y ) = \k x k y \ + E <k a - 2^/mf K -\ 

a—x 

(55) 

The 4>i a (x) is the wave function of the atoms at state p a 
on site i. 



Here G is the vector in the reciprocal lattice of the optical 
lattice, 4> s (k) is the Fourier transformation of the real 
space 4> s {r), which roughly behaves as 4> s {k) ~ const. 
at small k. The momentum distribution function of the 
superfluid phase of p— wave particles has been calculated 
inii. The result is complicated by the fact that in this 
phase there is large set of Z 2 degeneracies, as shown in 
(fM)) . The result of the calculation depends on which of 
the degenerate states the system is in. 

For a Mott Insulator phase of s— wave particles, the 
momentum distribution is 

n(k)/N~ \Mk)\ 2 - (58) 

For a Mott insulator of p— wave particles, one finds, 

n(k)/N ~ (a|^(fc)| 2 + b\<f> y (k)\ 2 ). (59) 

Here, the magnitudes of a and b depend on the relative 
abundance of p x and p y particles. At small momentum, 
4>x(k) ~ k x , (f> x (k) ~ k y . 

Although the algebraic bond liquid and the Mott In- 
sulator phases are qualitatively different, the momentum 
distribution functions are the same, 



n(k)/N ~ (\Mk)\ 2 n x + \^ y {k)\ 2 n v ). 



(60) 
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The reason for this is that the single particle correla- 
tion function is short-ranged in both the Mott insulator 
and the algebraic bond liquid nature, behaving in the 
latter case as (d^^d^) ~ Sij. Thus the momentum dis- 
tribution function cannot distinguish between these two 
insulating phases. 

C. Momentum density - momentum density 
correlator 

Recently it has been proposed that the momentum 
density - momentum density correlation function 

(n(R 1 )n(R 2 )) ~ {n Ql n Q2 ), (61) 

with Qi = mRi/t can be extracted in cold atom systems 
by measuring spatial noise correlations 21 . Remarkably, 
this technique has been successfully implemented to de- 
tect crystalline density order—. 

In the superfluid phase of p— wave particles, this cor- 
relator has been calculated inii, and the result again de- 
pends on which of the Z2 degenerate ground states the 
system is in. As we now show, this correlator can in prin- 
ciple be used to distinguish between the algebraic bond 
liquid and the Mott insulator phases. The calculation is 
easier to carry out for the normal ordered operator, 

{^l^M/N 2 = {n k n k ,)/N 2 + {n k )S(k - k')/N 2 , 

(62) 

and (n k )/N is given in (f60|) . The calculation of the cor- 
relation function in the algebraic bond liquid phase is 
lengthy but straight forward. At small k we find, 

(n k ng)/N 2 
~ a{h x k 2 x + n y k 2 )(n x k' 2 + n y k' 2 ) 
+ b{k x k' x n x + k y k' y n y ) 2 S(k -k' + G l ) 

i 

+c 1 (kl-k 2 y ){k' 2 -k' 2 )/N 
+a(k x k' x - k y k' y ) 2 F(k x - k' x , k y - k' y )/N 
-~fk 2 5(k-k')/N + 

(63) 

Here N is the total number of sites in the system, and 
the ellipses include all the terms smaller than 0(1/N). 
F is the function defined in (|55p. Vectors in the recip- 
rocal lattice of the optical lattice are denoted as G. The 
algebraic correlation function between flavor density give 
rise to the third and forth line in equation (|63[) . 

In the stripe ordered Mott Insulator phase defined in 
(|4ip , the momentum density - momentum density corre- 
lation function reads, 

(n k ng)/N 2 



~ a{n x k 2 x + n y k 2 y ){n x k' 2 + n y k' 2 ) 
+ ^ b(k x k' x n x + k y k' y n y ) 2 5{k - k' + G l ) 

i 

+ Y / ^(k X k' X ^k y k'y) 2 5(k-k' + G',) 

i 

-~fk 2 6(k -k')/N + •••, 

(64) 

where G" is the ordering wave vector of the stripe or- 
der. Therefore, this correlation function can be used to 
distinguish between the Mott insulator and the algebraic 
bond liquid phase. At the leading order, the two phases 
differ by the ordering wave vector G' . There are also 
qualitative differences in the 0{1/N) terms at small k. 

A more complicated correlation function which can in 
principle be used to directly measure the power law cor- 
relations in the algebraic bond liquid phase is (n\nf,,). 
We find, 

(n 2 n 2 ,)/N 4 ~ 0(1) + • • • + a(l/N 3 )k 2 x k 2 y k' 2 k' 2 
x {\k x -k' x \ 2 l^ K ^ + \k y ~k' y \ 2 /^ K ^ 

+ \k X + k'J 2 /^^ 1 + \ky + k'yl 2 ^^- 1 ). (65) 

The coefficient K is defined in the Gaussian action p4[) . 
and it directly encodes the power law decays in the alge- 
braic bond liquid phase. 

VIII. CONCLUSIONS AND EXTENSIONS 

In this work, a novel bosonic liquid phase has been 
predicted and explored in a system with strongly corre- 
lated bosons in different orbital levels of the wells of an 
optical lattice. In this algebraic bose liquid phase the 
correlation functions between two bond operators fall off 
as a power of the spatial separation. The reason for the 
existence of this new phase is the spatial anisotropy of 
the orbital wave functions: the p-level states extend and 
hop preferentially in one direction. Due to this orbital 
anisotropy, such a system in either d = 2, 3 can behave 
similarly to a one dimensional system. This effect has 
been called "dimensional reduction" 24 i 25 . From the the- 
oretical point of view, the phases studied in this work 
are new bosonic liquid phases in dimensions higher than 
one, which broadens our understanding of gapless "in- 
sulating" bosonic phases. Experimentally, it is an open 
challenge to see if such algebraic phases can be achieved 
in the cold atoms context. 

A possible extension of this work is to consider a 
fermionic version of the model studied herein. Fermionic 
atoms such as 40 K and 6 Li trapped in an optical lattice 
have been studied extensively recently, both theoretically 
and experimentally. Various behaviors have been pre- 
dicted - and some observed - in such fermionic atomic 
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systems. For instance, the normal fermi surface—, the 
condensate of molecules 40 K 2 and 6 Li 2 23 i 27 , the fermionic 
superfluid state with and without imbalanced fermion 
populations 2 ^ 2 ^ and the crossover between molecular 
condensate and the BCS like cooper pair condensate (for 
instance, see^). 

If fermionic atoms are pumped to p— level states, when 
interaction between atoms keeps the filling on each well 



close to n, we can again effectively study the low energy 
projected Hilbert space with the constraint. The ring 
exchange term similar to (fPT)|) is generated at low en- 
ergy. New fermionic phase different from all the previous 
studies is expected in this case, characterized by novel 
non-fermi liquid behaviors. 

The authors thank Eugene Demler for useful discus- 
sions. 
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